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ABSTRACT 

Observations and theory suggest that line driven winds from hot stars and luminous accretion disks adopt a 
unique, critical solution which corresponds to maximum mass loss rate. We analyze the numerical stability of the 
infinite family of shallow wind solutions, which resemble solar wind breezes, and their transition to the critical 
wind. Shallow solutions are sub-critical with respect to radiative (or Abbott) waves. These waves can propagate 
upstream through shallow winds at high speeds. If the waves are not accounted for in the Courant time step, 
numerical runaway results. The outer boundary condition is equally important for wind stability. Assuming pure 
outflow conditions, as is done in the literature, triggers runaway of shallow winds to the critical solution or to 
accretion flow. 
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1. Introduction 

Line driven winds (LDWs) occur in various astro- 
nomical objects, like OB and Wolf-Rayet stars, in ac- 
cretion disks in cataclysmic variables and, probably, in 
active galactic nuclei and luminous young stellar ob- 
jects. These winds are powered by absorption and re- 
emission of photosphcric continuum flux in numerous 
spectral transitions of C, N, O, Fc, etc. ions. 

Castor, Abbott, & Klein (1975; hereafter CAK) have 
analyzed the steady state Euler equation for LDWs. 
They found an infinite family of mathematical solu- 
tions, but only one, hereafter 'critical solution', which 
extends from the photosphere to arbitrary large radii. 
Other solutions either do not reach infinity or the pho- 
tosphere. The former solutions are called shallow and 
the latter ones steep. The unique, critical wind starts 
as the fastest shallow solution and switches smoothly 
to the slowest steep solution at the critical point. 

Observational support that LDWs adopt the criti- 
cal solution comes from measured terminal speeds (Ab- 
bott 1982). Furthermore, mass loss rates of supcrgiant 
winds are in general agreement with modified CAK the- 
ory (Lamers & Leitherer 1993; Puis et al. 1996). These 
measurements were recently extended to include galac- 
tic and extragalactic OB and A stars and central stars 
of planetary nebula (Kudritzki et al. 1999). 

Abbott (1980) put CAK theory in a complete anal- 
ogy to the solar wind and nozzle flows. The existence 
of a sonic point defines the unique, transsonic solutions 
for these flows, whereas the existence of a critical point 
for Abbott waves defines the unique, CAK solution 
for LDWs. Only from below this critical point, Ab- 
bott waves can propagate upstream towards the pho- 
tosphere. Above the critical point, they are advected 
outwards. Because Abbott waves generally propagate 
highly supersonically, the critical point of LDWs lies at 
much higher speeds than the sonic point. 

Abbott's (1980) analysis was challenged by Owocki 
& Rybicki (1986), who derived the Green's function for 
a pure absorption LDW. The Green's function gives 
correct signal speeds in presence of hydrodynamic in- 
stabilities. The inward signal speed in a pure absorp- 
tion line wind is the sound speed, and not the much 
larger Abbott speed, because photons propagate out- 
wards only. Owocki & Rybicki (1986) showed that a 
fiducial upstream signal, which still propagates inward 
at Abbott speed, must be interpreted as purely local 
Taylor series reconstruction. For a flow driven by scat- 
tering lines, however, Owocki & Puis (1999) find phys- 
ically relevant Abbott waves for a numerical Green's 
function. 

In the present paper, we further analyze the prop- 
erties of Abbott waves. We show that they are crucial 



for our understanding of stability of LDWs and must 
be included in the Courant time step. So far, time- 
dependent numerical simulations of LDWs from stars 
and accretion disks have ignored the ability of Abbott 
waves to communicate in the supersonic regime, which 
results in a numerical runaway. In particular, this run- 
away can lift the wind to the critical solution. 

The critical solution is also enforced by applying 
pure outflow boundary conditions. It is often argued 
that outflow boundary conditions are appropriate since 
LDWs are highly supersonic. Instead, they have to be 
super- abbottic. We show that shallow wind solutions, 
which correspond to solar wind breezes, are everywhere 
sub-abbottic. Hence, these solutions arc numerically 
destabilized by applying outflow boundary conditions. 

We formulate boundary conditions which render 
shallow solutions numerically stable. Those include 
non-reflecting Ricmann conditions for Abbott waves. 
By allowing for kinks in the velocity law, shallow solu- 
tions can be made globally admissible. 

2. The wind model 

In the CAK model for LDWs, both gravity and line 
force scale with r~ 2 . If the sound speed and hence the 
pressure forces are set to zero, this leads to a degen- 
eracy of the critical point condition, which is satisfied 
formally at every radius (Poe, Owocki, & Castor 1990). 
Thus, for this case, Abbott waves cannot propagate in- 
wards from any location in the wind. For finite sound 
speed, they creep inwards at small speed. Inclusion of 
the finite disk correction factor is much more relevant 
for LDWs than inclusion of pressure forces. With the 
finite disk included, the inward speed of Abbott waves 
below the critical point is significantly larger than the 
wind speed. Unfortunately, the finite disk correction 
factor depends on the (unknown) velocity law of the 
wind, which prevents a simple analysis of the wind dy- 
namics. 

We consider, therefore, a wind model which is ana- 
lytically feasible and yet prevents the (near-)degeneracy 
of the CAK point-star wind. (Especially, the latter 
leads to poor convergence of time-dependent numeri- 
cal schemes.) As a prototype, a vertical LDW from 
an isothermal, geometrically thin, non-self-gravitating 
accretion disk is assumed. The sound speed is set to 
zero. Keplerian rotation is assumed within the disk 
and angular momentum conservation above the disk. 
This reduces the flow problem to a 1-D, planar one. 
The radiative flux above an isothermal disk is roughly 
constant at small heights. On the other hand, the ver- 
tical gravity component along the wind cylinder is zero 
in the disk midplanc, grows linearly with z if z <C R 
(with R the footpoint radius in the disk), reaches a 
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maximum, and drops off at large z. To model the 
launch region of the wind and the gravity maximum, we 
choose g(z) — z/(l + z 2 ), with normalization GM = 1 
and R = 1, G being gravitational constant, and M is 
the mass of the central object. The different spatial 
dependence of flux and gravity results in a well-defined 
critical point in the flow. 

For constant radiative flux, the CAK line force be- 
comes g\ = C(p~ 1 dv/dz) a , where p and v are the den- 
sity and velocity, and a and C are constants. We choose 
a = 1/2. Typical values of a from NLTE calculations 
(Pauldrach 1987) and Kramers' law (Puis, Springmann, 
& Lennon 2000) are a = 0.4 to 0.7. The constant C can 
be expressed in terms of the unique, critical CAK mass 
loss rate. This is the maximum allowed mass loss rate 
for a stationary wind. We introduce new fluid variables 
m and w, 



where a subsript 'c' refers to the critical CAK solu- 
tion. For stationary, plane parallel winds, the continu- 
ity equation becomes pv = const, and m is a constant. 
In the following, the quantity w' = dw/dz = vdv/dz 
will play a central role. The velocity law, v(z), is ob- 
tained by quadrature of the wind 'solution' w'(z). For 
a = 1/2, the stationary Euler equation becomes, 

E = w' + g(z) - 2^g c w'/m = 0, (2) 

where we multiplied the nominator and denominator 
in the line force by v, and introduced g c = g(z c ), the 
gravity at the location of the CAK critical point. The 
dependence of w' on z is supressed in the equation. The 
constant C was determined as follows. Equation (||) is a 
quadratic equation in \fvJ . Consider the y/v/E plane. 
For arbitrary z and sufficiently small m, the crossings 
E = of the parabola E(y/v7) with the abscissa de- 
termine two wind solutions w'{z), termed shallow and 
steep. When m increases, the crossings approach each 
other, and merge at some maximum m > 1, beyond 
which no further solutions exist. Hence, the singularity 
condition dE/dw' = holds when two solutions merge, 
as is the case at the critical point. Together with E = 0, 
this gives w' c = g c and C — 2^Jg c p c v c . 

We add a comment which is important for the follow- 
ing. At every location z / z c , m > 1 for the parabola 
which meets the abscissa in just one point. The min- 
imum, m = 1, is reached at the unique critical point, 
z c , which is the How nozzle and determines the max- 
imum mass loss which can be driven from the photo- 
sphere to infinity. Solutions (actually, pairs of shallow 
and steep solutions) with m > 1 are called overloaded. 
They correspond to a choked Laval nozzle flow. These 



solutions become imaginary in a neigborhood of the 
critical point. We have shown elsewhere (Feldmeier 
& Shlosman 2000) how overloaded LDW solutions can 
be made globally defined, by letting them decelerate, 
dv/dz < 0, around the critical point. 

The position of the critical point is determined by 
the regularity condition, which is derived next. Since 
E = everywhere, dE/dz = holds too. Since 
dE/dw' c = 0, it follows that dE/dz c = 0, if |<| < oo. 
This latter condition singles out the unique critical 
point, since overloaded winds have kinks at the merging 
point of two solutions, \w"\ = oo. Hence, the critical 
point lies at the gravity maximum, dg/dz c = 0. For 
g(z) = z/(l + z 2 ), z c = 1 and g c = 1/2. The line force 
in (|[) becomes y/2w' /m. This holds for any solution, 
whether shallow or steep. Solving (|J) for Vu7, 

q(z) = y/2mw'{z) = 1 ± y/l - 2mg(z), (3) 

where we introduced a new variable q. The signs refer 
to steep (+) and shallow (— ) solutions. For m = 1, the 
square root in (||) vanishes at the critical point. For 
m > 1, the root becomes imaginary in a neighborhood 
of the critical point. The solutions w'(z) of (|3j) have a 
saddle topology in the zw' plane when m is varied from 
m < 1 to m > 1 . The critical point z c = 1 , w' c = 1 /2 is 
the saddle point. A variable radiative flux above a non- 
isothermal accretion disk leads instead to two different 
z c , above and below the gravity maximum (Feldmeier 
& Shlosman 1999). The lower is a saddle, the upper an 
extremum. 

Note that the CAK critical point of LDWs is a saddle 
in the z — vv' plane, whereas the sonic point of the solar 
wind is a saddle in the zv plane. In all other respects, 
we find a deep analogy between LDWs and the solar 
wind. Only the notion of a critical speed is replaced 
by a critical acceleration. Shallow solutions are the 
analogs to solar wind breezes. 

According to (||), shallow and steep solutions with 
m < 1 extend from z = to oo for the present model 
with zero sound speed. CAK found that steep solu- 
tions start supersonically at the wind base z — 0, which 
seems unphysical. In a spherically symmetric wind, 
shallow solutions do not extend to arbitrarily large z, 
because they cannot provide the required expansion 
work. These two results and the requirement for a 
continuous and differentiable solution implies that the 
wind has to pass through the critical point. It starts 
on the fastest shallow, and ends on the slowest steep 
solution. 

Figure 1 (which is adapted from Fig. 4 in Abbott 
1980 or Fig. 8.13 in Lamers & Cassinelli 1999) shows 
the LDW solution topology in the rv plane, for a fi- 
nite sound speed. In region (II), which is spatially the 
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most extended in real LDWs, shallow and steep solu- 
tions exist at each point (r, v). If z ^ z c , these solu- 
tion pairs may be overloaded. In the subsonic region 
(I), only shallow solutions (possibly overloaded) exist; 
steep solutions are necessarily supersonic. Correspond- 
ingly, only steep solutions (possibly overloaded) exist 
in region (III). No solutions exist in regions (IV) and 
(V). Note that the gap in overloaded solutions around 
z c can be bridged by an extended decelerating region, 
dv/dr < 0. 
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Fig. 1. - Solution topology of line driven winds, show- 
ing shallow, steep, and overloaded solutions. The bullet 
marks the CAK critical point. Roman numerals refer 
to cases (I) to (V) of CAK. Regions (I,II,III), and (V) 
are mutually separated by straight lines (sound speed 
v = a; Parker point r — rp); regions (II) and (IV) arc 
separated by the singular locus. 

The argument by CAK ruling out global, steep solu- 
tions is stringent, but the exclusion of shallow solutions 
seems too restrictive. The pressure force due to spheri- 
cal expansion becomes important only at very large dis- 
tances from the wind base. According to CAK, shallow 
solutions break down around 300 stellar radii for real- 
istic O star wind parameters. Most shallow solutions 
have by then overcome the local escape speed. Shal- 
low solutions can be globally defined by a slight gen- 
eralization of the CAK approach, allowing for negative 
v' = dv/dz in the line force <~ Vv*. This can be done 
by replacing v' by \v'\ or by max(u',0). The former 
expresses the 'blindness' of purely local line driving to 
flow acceleration or deceleration. The latter accounts 
for the nonlocal effect of shadowing in non-monotonic 
velocity laws, creating a resonance coupling. All radi- 
ation is assumed to be absorbed at the first resonance 
location. Realistically, the value of the line force lies 
between these extremes. The decelerating wind solu- 
tion, w' < 0, is unique, and does not split into shallow 



and steep branches. 

Hyperbolic differential equations allow for weak dis- 
continuities, which are represented in the zw' plane 
by jumps between different solution branches. Hence, 
a shallow solution may jump onto the decelerating 
branch. Even an overloaded solution can avoid becom- 
ing imaginary by jumping onto the decelerating branch 
in the vicinity of the critical point. This opens the 
possibility that global solutions can be composed in a 
piecewise manner. 

Jumps in w' introduce kinks in v{z). Such kinks have 
been found on various occasions in time-dependent 
LDW simulations. They may play a central role in 
the observed discrete absorption components in non- 
saturated P Cygni line profiles (Cranmer & Owocki 
1996). In LDWs from accretion disks in cataclysmic 
variables, a jump onto the decelerating branch even 
occurs for the critical solution, which would elsewhere 
not extend to infinity (Feldmeier & Shlosman 1999). 
Hence, the smoothness argument of CAK cannot be 
used to disqualify shallow solutions. 

3. Abbott waves 

3.1. Linear dispersion analysis 

The time-dependent continuity and Euler equation 
for the present wind model are, 



dp d(pv) _ 

dt + dz ~ u ' 

dv dv z 
dt dz 1 + z 2 



dv/dz 



(4) 
(5) 



In the Euler equation, v' > is assumed. The case v' < 
will be discussed later. To derive Abbott waves, small 
perturbations are applied, v = vq + v\ exp[i(fcz — tot)] 
and p = po+pi exp[i(fcz— cut)]. The subscript refers to 
a stationary solution, subscript 1 to the perturbation. 
Linearizing gives 
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(7) 



with non-dimensional \ = kv /v' and q — y / 2mow' . 
Setting the determinant of the system to zero gives the 
dispersion relation uj(k). In WKB approximation, \ ^ 
1, one finds for the phase speed, A , and growth rate 
of the inward (— ) mode, 



(8) 



A Q = Rc(w_)/fc = v ^1 - — ) , 
Im(w_) = 0, 
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in the observers frame, and for the outward (+) mode, 

A = Re(co + )/k = v , (9) 
Im(w+) = -e, 

with — e a small damping term of no further conse- 
quence. Along a shallow solution, qo < 1, and along 
a steep solution, go > 1. For shallow solutions with 
small mo and w' , qo is small, and the inward mode can 
propagate at arbitrarily large phase speed, from every 
location z. This result is new, since so far Abbott waves 
along shallow solutions were not considered ("a LDW 
has no analog to a solar breeze," Abbott 1980). Since 
for the CAK point star model every point is critical, 
qo = 1, and the inward mode is stagnating everywhere. 
Stability of non-WKB waves, \ = 0(1), is found from 
a numerical solution. 

Integrating (^|), the velocity of the critical solution 
at z c = 1 (the 'critical speed') is v w 0.61. On the other 
hand, a shallow solution with mo = 0.8 has v > 0.61 
for z > 1.51. This is no contradiction to the fact that 
the shallow solution is sub-abbottic. 

3.2. Characteristic analysis 

The above analysis holds for linear waves, and the 
wave speed is determined by the underlying stationary 
solution. General, non-linear waves can be derived from 
a characteristic analysis. The characteristic form of 
the above equations is, without further approximations 
(Feldmeier & Shlosman 2000), 



d_ d_ 

dt dz 

— A — 

dt dz 



dv 
dz 



ldv 
p dz 



ldg 
p dz ' 



(10) 
(11) 



where the Abbott speed, A, in the observers frame is 



A 



1 

1 - - 

q 



PcVc 



2p dv/dz' 



(12) 



The interpretation of (11) is unconventional, in that 
v' should now be considered as fundamental hydrody- 
namic field, instead of v. This happens because of the 
nonlinear dependence of E on v' (Courant & Hilbert 
1968). 

Hence, — pv' in the continuity equation and —g'/p in 
the Eulcr equation arc inhomogenous terms, since they 
do not contain derivatives of the hydrodynamic fields 
p and v'. The advection operators are d/dt + vd/dz 
and d/dt + Ad/dz, with characteristic speeds v and A. 
The former should be read as v + a in zero sound speed 
approximation. The Riemann invariants p and v'/p 
correspond to wave amplitudes. In a frame moving at 



v, the amplitude p is constant, except for sources and 
sinks —pv'. In an isothermal gas, p ~ p, and the wave 
amplitude is indeed that of a pressure wave. In a frame 
moving at A, the amplitude v'/p is constant, except for 
sources and sinks —g'j p. The Sobolev optical depth is 
~ (v'/p) , indicating that the (— ) mode is a radiative 
wave. 

Even for waves of infinitesimal amplitude, v' may 
deviate strongly from v' , if the wavelength A is short. 
Since A ~ (t/) -1 / 2 , short-scale Abbott waves are highly 
dispersive. For sufficiently small A, v' becomes nega- 
tive at certain wave phases. The characteristic equa- 
tions (10, 11) hold also for v' < 0, but A has different 
meaning then. For a line force ~ the general A 

is 



A 



1 

it — 

5± 



(13) 



where the lower sign applies for v' < 0, and in g_ 
one has to use —v'. Trivially, Abbott waves are ad- 
vected outwards in super-abbottic flow. But that Ab- 
bott waves propagate downstream as a (+) mode in re- 
gions where v' < is probably their single, most unique 
property. 

For a line force ~ ^/max(0, v') then, A = v if v' < 0. 
This is because the line force drops out of the Eu- 
ler equation, and both wave modes become ordinary 
sound. (Actually A — v — a then, for upstream propa- 
gating sound.) In this case one has along a character- 
istic with coordinate s, 
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P(s) = v'( s ) = 

p(0) v'(Q) l + v'(0)s' 



(14) 



as WKB wave solution of (10, 11). This describes gas 
dynamical simple waves. Since v'(0) < 0, they steepen 
and break (v' = -co). The derivation of simple Abbott 
waves is beyond the scope of this paper. 

4. Boundary conditions for line driven winds 

We discuss now the outer boundary conditions for 
shallow solutions. Abbott waves can enter through the 
outer boundary, hence, pure outflow boundary condi- 
tions are not appropriate for shallow solutions. 

As hydrocode we use a standard Eulerian scheme on 
a staggered mesh. Advection is solved in the conser- 
vative control volume approach, using monotonic van 
Leer (1977) interpolation. The sound speed is set ex- 
actly to zero, and a small amount of artificial viscosity 
is added to handle occasional shock fronts. Details on 
the scheme can be found in Reile & Gehren (1991) and 
Feldmeier (1995). 

In a first step, we consider simple boundary condi- 
tions for p and v which still preserve essential features 
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of Abbott waves. These boundary conditions render 
shallow solutions stable, and allow small perturbations 
to leave the grid. They are, 

p(0) = const, v(0)=v(l), 
p(I) = p (I - 1), v(I) = const, (15) 

where and I are the mesh indices of the inner and 
outer boundary, respectively. The conditions (15) are 
motivated as follows. (The reader may jump to the end 
of this paragraph.) On a characteristic curve [z(s),t(s)] 
which leaves the mesh through a boundary, the charac- 
teristic equation is solved on the boundary via extrap- 
olation, using one-sided, interior differentials. Such an 
upwinded scheme is stable (Steger & Warming 1981). 
For shallow winds, the p wave leaves through the outer 
boundary if v > 0, hence p is extrapolated from the 
interior in (15). Zeroth order extrapolation is used 
for stability reasons. Abbott waves enter the mesh 
through the outer boundary, hence a boundary con- 
dition must be applied. We choose v(I) = const. (To 
prevent standing waves, v(I) — m /p(I) is often more 
appropriate; but is also more susceptible to runaway 
than v(I) = const.) At the inner boundary, the argu- 
ment proceeds correspondingly, with interchanged roles 
of the waves. We find that using p(0) — p(X) and 
v(0) = const instead destabilizes shallow solutions. 

4.1. The Courant time step 

In time-dependent hydrodynamic simulations pub- 
lished so far, it is customary to insert the sound speed 
in the Courant time step, 

At = a mm (16) 
\v ± a\ 

with Courant number a < 1, and the minimum has 
to be taken over the mesh. We see now that Abbott 
waves, as characteristic inward wave mode, have to be 
included in the Courant time step, 

( Az Az \ 
At = a mm -, , 17) 

\\v - A\' v + a J 

where A is in the fluid frame now. Note that A changes 
sign with v'. For sufficiently small v' > and for ar- 
bitrary v' < 0, Abbott waves determine the time step. 
At velocity plateaus where the wind velocity is more 
or less constant, A — > — oo from (11), and the Courant 
time step — * 0. This is an artifact of Sobolev approx- 
imation in which the line radiation force tends to zero 
in the absence of velocity gradients. In practice, in 
calculating the Courant time step we do not allow q 
to drop below a certain minimum value, usually 10 . 
This corresponds to a « 10 -4 for an ordinary Courant 



time step not including Abbott waves. Our results are 
not sensitive to the value of this minimum q, if the lat- 
ter is sufficiently small. Furthermore, velocity plateaus 
quickly acquire a tilt and propagate at finite speed. 

4.2. Non-reflecting boundary conditions 

We can then formulate non-reflecting boundary con- 
ditions for the Riemann invariants p and v 1 / p, which 
annihilate incoming waves. This prevents boundary re- 
flection of waves which originate on the interior mesh. 
To annihilate a wave, its amplitude (the Riemann in- 
variant) is kept constant in time at the boundary, 
d(v'/p)/dt — and dp/dt = for Abbott and sound 
waves, respectively (Hedstrom 1979; Thompson 1987). 
One may say that non-reflecting boundary conditions 
drive the numerical solution towards a neighboring, sta- 
tionary solution without waves. The technical details of 
implementing these boundary conditions are discussed 
in Appendix A. 

5. Stability of solutions 

We discuss here the numerical stability of shal- 
low and steep solutions, depending on the appropriate 
Courant time step and boundary conditions. 

5.1. Stability of steep solutions 

We start with steep solutions, and show that they 
are unconditionally stable, even when simplified bound- 
ary conditions analogous to (15) are used. According 
to (§), q > 1 everywhere for steep solutions, and Ab- 
bott waves can propagate only towards larger z. At 
the outer boundary, extrapolation p(I) = p(I — 1) and 
v(I) = v(I — 1) is therefore appropriate. At the inner 
boundary, two conditions have to be specified. The 
obvious choices are either p — const, v — const or 
p = const, v' — const. The first set fixes the mass 
loss rate, the second establishes non-reflecting bound- 
ary conditions: the wave amplitudes are p and v' / p, and 
the latter boundary conditions keep them constant in 
time, meaning there are no waves. 

We find that steep solutions are stable for both types 
of boundary conditions. Even when the initial condi- 
tions differ profoundly from a steep wind, the numeri- 
cal code converges to the latter. Furthermore, we per- 
formed tests with an explicit, harmonic perturbation 
at fixed Eulerian position in the wind. Even if the per- 
turbation amplitude reached 100%, the wind remained 
on average on a steep solution, which can, therefore, be 
considered as unconditionally stable. 
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5.2. Stability and runaway of shallow solutions 

A physically more relevant question is for the sta- 
bility of shallow solutions. In a first step, we use an 
analytic, shallow wind as initial conditions. Table [l] 
shows the runaway of this shallow to the critical ve- 
locity law. Abbott waves are not accounted for in the 
Courant time step, and pure outflow boundary condi- 
tions p(I) = p(I — 1) and v(I) = v(I — 1) are applied. 
This corresponds to the typical procedure adopted in 
the literature. The runaway starts at the outer bound- 
ary and generates inward propagating Abbott waves. 
The velocity law evolves towards larger speeds, until 
the critical solution is reached. 

The mechanism of the runaway seems as follows. 
The Courant condition for Abbott waves is violated 
in the outer wind, where the velocity law is flat and 
the Abbott speed is large. This causes numerical run- 
away. Opposed to standard gas dynamics, the LDW 
runaway does not crash the simulation. We speculate 
that exponential growth in v creates large \v'\, which 
imply small Abbott speeds. The Courant condition is 
no longer violated, and runaway stops. Abbott waves 
are also excited at the outer boundary, by switching 
from the interior to the boundary scheme. Upstream 
propagating Abbott waves are inconsistent with the as- 
sumed outflow boundary conditions. The waves propa- 
gate to the wind base and drive the inner wind towards 
the critical solution, which is indeed consistent with 
outflow boundary conditions. 

Note in Table that Abbott waves can propagate 
inwards from the outer boundary, zm = 4, until the 
critical velocity law is reached (indicated by bold num- 
bers in the table). This is no contradiction to the crit- 
ical point being located at z c = 1: the outer, shallow 
portion of the velocity law (the numbers in italics) is 
sub-abbottic, even when the inner, steep velocity law 
is already super-abbottic. 

To prevent numerical runaway of shallow solutions, 
we apply non-reflecting Riemann boundary conditions 



Table 1: Numerical runaway of a shallow velocity law. 
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0.8 


1.2 


1.4 


1.4 


1.4 


2.5 


0.8 


0.9 


1.3 


1.5 


1.8 


1.8 


3.0 


0.8 


1.1 


1.3 


1.5 


2.1 


2.1 


3.5 


0.8 


1.1 


14 


1.6 


2.2 


2.4 


4.0 


0.8 


1.2 


1.4 


1.6 


2.3 


2.7 



and include Abbott waves in the Courant time step. 
As initial conditions, an arbitrary (for example, linear) 
velocity law is used. Abbott waves are again exited at 
the outer boundary and propagate inwards. However, 
the simulation relaxes quickly to a shallow solution. 
The m-value depends on the initial data. Even for an 
initial model with height-dependent mass flux, we find 
quick convergence to a shallow solution. 

Further details on the stability of shallow solutions 
are given in Appendix B. To summarize this section, 
both improper outer boundary conditions and Courant 
conditions which do not account for Abbott waves are 
responsible for numerical runaway. 

6. Summary 

We present a simplified model for line driven winds 
from stars and accretion disks which avoids the r~ 2 
degeneracy of the CAK model. Radiative or Abbott 
waves are derived, both from a dispersion and a char- 
acteristic analysis, and for all possible solutions to the 
Euler equation, i.e., shallow, critical, steep and over- 
loaded ones. 

We find that Abbott waves can propagate upstream, 
towards the photosphere, from any position along a 
shallow wind solution. Hence, for the latter, they can 
enter the calculational grid through the outer bound- 
ary The wave propagation speed depends on the ve- 
locity slope of the wind. Abbott's (1980) result that 
these waves are creeping, i.e., have only slighly nega- 
tive inward speeds, holds globally only for the almost 
degenerate CAK point-star model. In more realistic 
wind models, Abbott waves can limit the Courant time 
step. 

Both the neglect of Abbott waves in the Courant 
time step and in the boundary conditions leads to nu- 
merical runaway towards the critical wind solution of 
maximum mass loss rate or to accretion flow. If, in- 
stead, incoming waves are annihilated at the outer 
boundary and the correct Courant timestep is used, 
shallow solutions are stable. 

It is a pleasure to thank Janet Drew, Leon Lucy, 
Stan Owocki, and Joachim Puis for frequent discus- 
sions. This work was supported in part by NASA 
grants NAG 5-10823, NAG5-3841, WKU-522762-98-6 
and HST GO-08123.01-97A to I.S., which are grate- 
fully acknowledged. 

A. Non-reflecting boundary conditions 

We adopt the following procedure to calculate p and 
v on a boundary, (i) If the p wave leaves the grid, 
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p = —(pv)' is calculated using one-sided, interior dif- 
ferentials, (ii) If the p wave enters the grid, p = is set 
(wave annihilation; see the main text), (iii) If the v' j p 
wave leaves the grid, v'/v' is calculated from the Euler 
equation (11) in characteristic form, using one-sided, 
interior differentials. Note that v" appears here, (iv) If 
v'/p enters the grid, v'/v' = p/p is set, assuming WKB 
waves, (v) The new values for p and v on the boundary 
are found by integrating p and v' using a time-explicit 
scheme. 

In the latter step (v), one actually needs u, not v' . 
To integrate from time step n — 1 to n at the outer 
boundary, using v' , we write 



v'At. 



(Al) 



Using v' n — (vf — u"_ 1 )/Az, and similarly for v' n , 
this becomes, 



»p = v 11 - 1 + - v'izl + v'AzAt, 



(A2) 



whence v T n depends on v™~Z\ ■ Due to frequent variable 
updating in an operator splitting scheme during the 
time step, u^-i wou ld have to be stored as an extra 
variable. This seems not desirable, and causes indeed 
numerical runaway of sha llow solutions. Fortunately, a 
simple approximation to ( |A2| ) gives satisfactory results, 
namely setting v"Z± = v 7 }_ l , or 



v'AzAt. 



(A3) 



B. Numerical stability of shallow solutions 

For numerical stability tests, we assume the follow- 
ing parameters, if not otherwise stated: 200 equidis- 
tant grid points from z = 0.1 to z = 4. Below 
z = 0.1, the stationary wind speed is very small. Since 
5v/Aq = Sp/po, one has vq <C Aq for an inner bound- 
ary z m — > 0. Even small perturbations Sp cause then 
negative speeds at z m . This alters the direction of the 
p characteristic, and often causes numerical problems. 
A shallow solution with to = 0.8 serves as start model. 
The sound speed is set to zero. The flow time from 
z = 0.1 to z = 4 is 12.4 for a shallow solution with 
to = 0.8, and 9.1 for the critical solution (in the units 
specified in the main text). For to = 0.8, linear Abbott 
waves propagate from z = 4 to z = 0.1 in a time 6.0. 
All simulations are evolved to time 50. Table B2 lists 
results from simulations using different outer bound- 
ary conditions and Courant time steps. On the inner 
boundary, conditions (15) were used. 

Consider first the line force ~ \/|t>'|. The table 



time is mandatory to prevent runaway (runs 1-5; 'shal- 
low' in the column for the outer boundary conditions 
refers to (15). For small Courant numbers, Abbott 
waves are, to a degree, already accounted for by an or- 
dinary time step. Runaway may then be prevented, as 
in run 2, though the solution does oscillate at all times. 
Shifting the outer boundary outwards, the runaway oc- 
curs again (run 5). If Abbott waves are accounted for 
in the time step, outflow boundary conditions (extrap- 
olation of p and v) do not necessarily lead to runaway. 
Instead, a stable, shallow solution may be maintained 
(runs 6, 9; also 2). However, the solution oscillates 
at all times. Shallow solutions become more unstable 
when the outer boundary is shifted outwards, into the 
shallow part of the velocity curve where Abbott waves 
are easily excited. Finally, the control run 10 shows 
that the initial shallow solution is numerically stable if 
Abbott waves are included in the time step and bound- 
ary conditions (15) are used. 

For a line force ~ yj max(0, v'), the scheme is even 
more susceptible to numerical runaway. All cases with 
outflow boundary conditions undergo runaway. How- 
ever, for the boundary condition v(I) = v(I — 1), the 
solution no longer tends towards the critical solution. 
Instead v(I) drops to zero, and with it the interior wind 
velocity. The outcome of the simulation is not clear, as 
we have not formulated appropriate line driven accre- 
tion boundary conditions. The reason for the drop in 
v(I) is that an accidental v' < at the outer bound- 
ary implies zero line force. This causes further veloc- 
ity drop and v' becomes more negative. Abbott waves 
carry this wind breakdown to the interior mesh. It is 
not clear whether this type of line force should actually 
be applied in the proximity of the outer boundary. 
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Table B2 

Runaway in simulations applying different boundary conditions and 

courant time steps. 



run 




grid 


line 


Abbott 


Courant 


boundary 


result 






points 


force 


time 


number 


condition 




1 


4 


200 


abs(u') 


no 


0.5 


outflow 


runaway 


2 


4 


200 


abs(u') 


no 


0.1 


outflow 


m = 0.7; oscill. 


3 


4 


200 


abs(u') 


no 


0.5 


shallow 


runaway 


4 


4 


200 


abs(u') 


no 


0.1 


shallow 


runaway 


5 


10 


200 


abs(u') 


no 


0.1 


outflow 


runaway 


6 


4 


200 


abs(v') 


yes 


0.5 


outflow 


m = 0.5; oscill. 


7 


10 


200 


abs(u') 


yes 


0.5 


outflow 


slow runaway 


8 


10 


500 


abs(V) 


yes 


0.5 


outflow 


slow runaway 


9 


4 


80 


abs(V) 


yes 


0.5 


outflow 


m = 0.5; oscill. 


10 


4 


200 


abs(w') 


yes 


0.5 


shallow 


stable 


11 


4 


200 


max(0, v') 


no 


0.5 


outflow 


accretion 


12 


4 


200 


max(0, v') 


no 


0.1 


outflow 


accretion 


13 


4 


200 


max(0, v') 


no 


0.5 


shallow 


stable, oscill. 


14 


4 


200 


max(0, v') 


no 


0.1 


shallow 


stable, oscill. 


15 


10 


200 


max(0, v') 


no 


0.1 


outflow 


accretion 


16 


4 


200 


max(0, v') 


yes 


0.5 


outflow 


accretion 


17 


10 


200 


max(0, v') 


yes 


0.5 


outflow 


accretion 


18 


10 


500 


max(0, v') 


yes 


0.5 


outflow 


accretion 


19 


4 


80 


max(0, v') 


yes 


0.5 


outflow 


accretion 


20 


4 


200 


max(0, v') 


yes 


0.5 


shallow 


stable 
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